﻿// SyntheticDataGenerator.Generators.UniformGenerator
//
// (c) 2011 Arthur Pitman
//
// Part of Market-Basket Synthetic Data Generator
// A C# adaptation of the IBM Quest Market-Basket Synthetic Data Generator
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// (Version 2) as published by the Free Software Foundation.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public
// License (Version 2) along with this program; if not, write to 
// the Free Software Foundation, Inc., 
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
// or see <http://www.gnu.org/licenses/>.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SyntheticDataGenerator.Generators
{
	/// <summary>
	/// Generates uniformly distributed random numbers
	/// 
	/// See Park, Stephen K. and Keith W. Miller. 
	/// “Random Number Generators: Good Ones Are Hard to Find”. 
	/// Communications of the ACM, October 1988, p. 1192.
	/// 
	/// See "Numerical Recipes in C: The Art of Scientific Computing" by Press et al.
	/// </summary>
	public class UniformGenerator : IGenerator
	{
		protected const int IA = 16807;
		protected const int IM = 2147483647;
		protected const double AM = 1.0 / IM;
		protected const int IQ = 127773;
		protected const int IR = 2836;
		protected const int MASK = 123459876;

		protected const int NTAB = 32;
		protected const int NDIV = (1 + (IM - 1) / NTAB);
		protected const double EPS = 1.2e-7;
		protected const double RNMX = (1.0 - EPS);

		protected int idum;
		protected int iy;
		protected int[] iv = new int[NTAB];

		/// <summary>
		/// Initializes a new <see cref="UniformGenerator"/> instance with a random seed
		/// </summary>
		public UniformGenerator()
		{
			iy = 0;
			idum = SeedGenerator.Next();
		}

		/// <summary>
		/// Initializes a new <see cref="UniformGenerator"/> instance with specific seed
		/// </summary>
		/// <param name="seed"></param>
		public UniformGenerator(int seed)
		{
			iy = 0;
			idum = seed;
		}

		/// <summary>
		/// Initializes the internal state of the <see cref="UniformGenerator"/>
		/// </summary>
		protected void Initialize()
		{
			// Be sure to prevent idum = 0
			if (-idum < 1)		
				idum = 1;
			else
				idum = -idum;

			// load the shuffle table (after 8 warmups)
			for (int index = NTAB + 7; index >= 0; index--)
			{
				int x = idum / IQ;
				idum = IA * (idum - x * IQ) - IR * x;

				// fix overflows by Schrage's method
				if (idum < 0)
					idum += IM;

				if (index < NTAB)
					iv[index] = idum;
			}
			iy = iv[0];
		}

		/// <summary>
		/// Generates a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint values)
		/// Random number generator based on that of Park and Miller
		/// Includes Bays-Durham shuffle and added safeguards
		/// Set idum negative to initialize; after this do notalter idum between successive deviates
		/// RNMX approximates the largest floating value less than 1.
		/// </summary>
		/// <returns>Uniform random deviate</returns>
		public double Next()
		{       
			// Check if initialization should be performed
			if ((idum <= 0) || (iy == 0))
				Initialize();

			// Start here when not initializing.
			int x = idum / IQ;			  

			// Compute idum = (IA*idum) % IM without
			idum = IA * (idum - x * IQ) - IR * x; 

			// Fix overflows by Schrage's method
			if (idum < 0) 
				idum += IM;

			// Calculate index in the range 0..NTAB-1
			int index = iy / NDIV;			  

			// Output previously stored value and and refull the shuffle table.
			iy = iv[index];            
			iv[index] = idum;
						
			// Limit to endpoint values
			double value = iy * AM;
			if (value > RNMX) 	  
				value = RNMX;
			
			return value;
		}

		/// <summary>
		/// Returns k distinct uniformly distributed random integers between 0 and (n-1), sorted in ascending order
		/// </summary>
		/// <param name="n"></param>
		/// <param name="k"></param>
		/// <returns></returns>
		public IEnumerable<int> ChooseRandomSubset(int n, int k)
		{           
			KeyValuePair<int, double>[] items = new KeyValuePair<int, double>[n];
			for (int i = 0; i < n; i++)
			{
				items[i] = new KeyValuePair<int, double>(i, Next());
			}

			// randomize the order of the items
			var sortedItems = from x in items
							  orderby x.Value ascending
							  select x.Key;

			// take first k and resort them 
			var subset = from x in sortedItems.Take(k)
						 orderby x ascending
						 select x;

			return subset;
		}
	}
}
